#!/usr/bin/env python
# -*- coding: utf-8 -*-

import shared

# If we are presented with the first k terms of a sequence it is
# impossible to say with certainty the value of the next term, as
# there are infinitely many polynomial functions that can model the
# sequence.

# As an example, let us consider the sequence of cube numbers. This is
# defined by the generating function,
# u_n = n^3: 1, 8, 27, 64, 125, 216, ...

# Suppose we were only given the first two terms of this
# sequence. Working on the principle that "simple is best" we should
# assume a linear relationship and predict the next term to be 15
# (common difference 7). Even if we were presented with the first
# three terms, by the same principle of simplicity, a quadratic
# relationship should be assumed.

# We shall define OP(k, n) to be the nth term of the optimum
# polynomial generating function for the first k terms of a
# sequence. It should be clear that OP(k, n) will accurately generate
# the terms of the sequence for n <= k, and potentially the first
# incorrect term (FIT) will be OP(k, k+1); in which case we shall call
# it a bad OP (BOP).

# As a basis, if we were only given the first term of sequence, it
# would be most sensible to assume constancy; that is, for n >= 2,
# OP(1, n) = u_1.

# Hence we obtain the following OPs for the cubic sequence:
# OP(1, n) = 1 	1, 1, 1, 1, ...
# OP(2, n) = 7n-6 	1, 8, 15, ...
# OP(3, n) = 6n^2-11^n+6      	1, 8, 27, 58, ...
# OP(4, n) = n^3 	1, 8, 27, 64, 125, ...

# Clearly no BOPs exist for k > 4.

# By considering the sum of FITs generated by the BOPs (indicated in
# red above), we obtain 1 + 15 + 58 = 74.

# Consider the following tenth degree polynomial generating function:

# un = 1 - n + n^2 - n^3 + n^4 - n^5 + n^6 - n^7 + n^8 - n^9 + n^10

# Find the sum of FITs for the BOPs.

expected = 37076114526

def u(n):
    return sum(((-n)**p for p in range(11)))

def equations(series, k):
    eqs = []
    for x in range(1, k+1):
        eq = []
        for j in range(k):
            eq.append(x**j)
        eq.append(series[x-1])
        eqs.append(eq)
    return eqs
    

def solve_equations(equations):
    coefficients = [0] * (len(equations[0])-1)
    for equation_number in range(len(equations)-1):
        first_eq = equations[equation_number]
        for eq in range(equation_number+1, len(equations)):
            multiplier = float(equations[eq][equation_number])/first_eq[equation_number]
            for i in range(equation_number, len(first_eq)):
                equations[eq][i] -= multiplier * first_eq[i]
            #print equations[eq], multiplier
    # we should now have the coefficients in upper-triangular form - walk back and
    # start picking out the individual coefficients

    for eq in range(len(equations)-1, -1, -1):
        # clear out higher coeffs
        for i in range(eq, len(equations)):
            equations[eq][-1] -= coefficients[i] * equations[eq][i]

        # figure out this coefficient
        coefficients[eq] = float(equations[eq][-1]) / equations[eq][eq]
    return coefficients


def solve():

    total = 0
    
    biggest_k = 10
    series = [u(n) for n in range(1, biggest_k+1)]
    #series = [n**3 for n in range(1, biggest_k+1)]
    for k in range(1, biggest_k+1):
        matrix = equations(series, k)
        #print 'matrix', matrix
        coefficients = solve_equations(matrix)
        #print 'coefficients', coefficients
        fit = sum(((k+1) ** n * coefficients[n] for n in range(len(coefficients))))
        total += fit

    return int(total)
        
#     for i in range(1, 11):
#         print i, u(i)
#     pass

